home *** CD-ROM | disk | FTP | other *** search
/ Ham Radio 2000 #2 / Ham Radio 2000 - Volume 2.iso / HAMV2 / MISC / DTMFF110 / REALFFT.C < prev    next >
Encoding:
C/C++ Source or Header  |  1997-06-10  |  4.7 KB  |  180 lines
  1. /*
  2.  *     Program: REALFFT.C Author: Philip VanBaren  Date: 2 September 1993
  3.  *
  4.  * Description: These routines perform an FFT on real data.
  5.  *        On a 486/33 compiled using Borland C++ 3.1 with full
  6.  *        speed optimization and a small memory model, a 1024 point
  7.  *        FFT takes about 16ms.
  8.  *        This code is for integer data, but could be converted
  9.  *        to float or double simply by changing the data types
  10.  *        and getting rid of the bit-shifting necessary to prevent
  11.  *        overflow/underflow in fixed-point calculations.
  12.  *
  13.  *  Note: Output is BIT-REVERSED! so you must use the BitReversed to
  14.  *      get legible output, (i.e. Real_i = buffer[ BitReversed[i] ]
  15.  *                    Imag_i = buffer[ BitReversed[i]+1 ] )
  16.  *      Input is in normal order.
  17.  *
  18.  *  Copyright (C) 1995    Philip VanBaren
  19.  */
  20.  
  21. /*
  22.  * For Borland C++/TASM you may define the flag ASMFFTCODE in order to use the
  23.  * assembly code for the FFT.  This will skip the compilation of the C version
  24.  * of the routine; you must be sure to add realffta.asm to your makefile
  25.  */
  26.  
  27. #include <math.h>
  28. #include <stdlib.h>
  29. #include <stdio.h>
  30.  
  31. int          *BitReversed;
  32. short        *SinTable;
  33. int           Points = 0;
  34.  
  35. /*
  36.  *  Initialize the Sine table and Twiddle pointers (bit-reversed pointers)
  37.  *  for the FFT routine.
  38.  */
  39. void 
  40. InitializeFFT( int fftlen )
  41. {
  42.   int           i;
  43.   int           temp;
  44.   int           mask;
  45.  
  46.   /*
  47.    * FFT size is only half the number of data points The full FFT output can
  48.    * be reconstructed from this FFT's output. (This optimization can be made
  49.    * since the data is real.)
  50.    */
  51.   Points = fftlen;
  52.  
  53.   if ( ( SinTable = ( short * ) malloc( Points * sizeof( short ) ) ) == NULL )
  54.   {
  55.     puts( "Error allocating memory for Sine table." );
  56.     exit( 1 );
  57.   }
  58.   if ( ( BitReversed = ( int * ) malloc( Points / 2 * sizeof( int ) ) ) == NULL )
  59.   {
  60.     puts( "Error allocating memory for BitReversed." );
  61.     exit( 1 );
  62.   }
  63.  
  64.   for ( i = 0; i < Points / 2; i++ )
  65.   {
  66.     temp = 0;
  67.     for ( mask = Points / 4; mask > 0; mask >>= 1 )
  68.       temp = ( temp >> 1 ) + ( i & mask ? Points / 2 : 0 );
  69.  
  70.     BitReversed[i] = temp;
  71.   }
  72.  
  73.   for ( i = 0; i < Points / 2; i++ )
  74.   {
  75.     register double s, c;
  76.     s = floor( -32768.0 * sin( 2 * M_PI * i / ( Points ) ) + 0.5 );
  77.     c = floor( -32768.0 * cos( 2 * M_PI * i / ( Points ) ) + 0.5 );
  78.     if ( s > 32767.5 )
  79.       s = 32767;
  80.     if ( c > 32767.5 )
  81.       c = 32767;
  82.     SinTable[BitReversed[i]] = ( short ) s;
  83.     SinTable[BitReversed[i] + 1] = ( short ) c;
  84.   }
  85. }
  86.  
  87. /*
  88.  *  Free up the memory allotted for Sin table and Twiddle Pointers
  89.  */
  90. void 
  91. EndFFT(  )
  92. {
  93.   free( BitReversed );
  94.   free( SinTable );
  95.   Points = 0;
  96. }
  97.  
  98. /* Include this only if you don't use the REALFFTA.ASM routine */
  99. /* This is for DOS only */
  100.  
  101. #ifndef ASMFFTCODE
  102.  
  103. short        *A, *B;
  104. short        *sptr;
  105. short        *endptr1, *endptr2;
  106. int          *br1, *br2;
  107. long          HRplus, HRminus, HIplus, HIminus;
  108.  
  109. /*
  110.  *  Actual FFT routine.     Must call InitializeFFT(fftlen) first!
  111.  */
  112. void 
  113. RealFFT( short *buffer )
  114. {
  115.   int           ButterfliesPerGroup = Points / 4;
  116.  
  117.   endptr1 = buffer + Points;
  118.  
  119.   /*
  120.    * Butterfly: Ain-----Aout \ / / \ Bin-----Bout
  121.    */
  122.  
  123.   while ( ButterfliesPerGroup > 0 )
  124.   {
  125.     A = buffer;
  126.     B = buffer + ButterfliesPerGroup * 2;
  127.     sptr = SinTable;
  128.  
  129.     while ( A < endptr1 )
  130.     {
  131.       register short sin = *sptr;
  132.       register short cos = *( sptr + 1 );
  133.       endptr2 = B;
  134.       while ( A < endptr2 )
  135.       {
  136.     long          v1 = ( ( long ) *B * cos + ( long ) *( B + 1 ) * sin ) >> 15;
  137.     long          v2 = ( ( long ) *B * sin - ( long ) *( B + 1 ) * cos ) >> 15;
  138.     *B = ( *A + v1 ) >> 1;
  139.     *( A++ ) = *( B++ ) - v1;
  140.     *B = ( *A - v2 ) >> 1;
  141.     *( A++ ) = *( B++ ) + v2;
  142.       }
  143.       A = B;
  144.       B += ButterfliesPerGroup * 2;
  145.       sptr += 2;
  146.     }
  147.     ButterfliesPerGroup >>= 1;
  148.   }
  149.   /*
  150.    * Massage output to get the output for a real input sequence.
  151.    */
  152.   br1 = BitReversed + 1;
  153.   br2 = BitReversed + Points / 2 - 1;
  154.  
  155.   while ( br1 <= br2 )
  156.   {
  157.     register long temp1, temp2;
  158.     short         sin = SinTable[*br1];
  159.     short         cos = SinTable[*br1 + 1];
  160.     A = buffer + *br1;
  161.     B = buffer + *br2;
  162.     HRplus = ( HRminus = *A - *B ) + ( *B << 1 );
  163.     HIplus = ( HIminus = *( A + 1 ) - *( B + 1 ) ) + ( *( B + 1 ) << 1 );
  164.     temp1 = ( ( long ) sin * HRminus - ( long ) cos * HIplus ) >> 15;
  165.     temp2 = ( ( long ) cos * HRminus + ( long ) sin * HIplus ) >> 15;
  166.     *B = ( *A = ( HRplus + temp1 ) >> 1 ) - temp1;
  167.     *( B + 1 ) = ( *( A + 1 ) = ( HIminus + temp2 ) >> 1 ) - HIminus;
  168.  
  169.     br1++;
  170.     br2--;
  171.   }
  172.   /*
  173.    * Handle DC bin separately
  174.    */
  175.   buffer[0] += buffer[1];
  176.   buffer[1] = 0;
  177. }
  178.  
  179. #endif
  180.